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ABSTRACT 

Measurements of the redshift-space galaxy clustering have been a prolific source of 
cosmological information in recent years. Accurate covariance estimates are an essen¬ 
tial step for the validation of galaxy clustering models of the redshift-space two-point 
statistics. Usually, only a limited set of accurate N-body simulations is available. Thus, 
assessing the data covariance is not possible or only leads to a noisy estimate. Further, 
relying on simulated realisations of the survey data means that tests of the cosmol¬ 
ogy dependence of the covariance are expensive. With these points in mind, this work 
presents a simple theoretical model for the linear covariance of anisotropic galaxy 
clustering observations with synthetic catalogues. Considering the Legendre moments 
(‘multipoles’) of the two-point statistics and projections into wide bins of the line-of- 
sight parameter (‘clustering wedges’), we describe the modelling of the covariance for 
these anisotropic clustering measurements for galaxy samples with a trivial geometry 
in the case of a Gaussian approximation of the clustering likelihood. As main result of 
this paper, we give the explicit formulae for Fourier and configuration space covariance 
matrices. To validate our model, we create synthetic HOD galaxy catalogues by popu¬ 
lating the haloes of an ensemble of large-volume N-body simulations. Using linear and 
non-linear input power spectra, we find very good agreement between the model pre¬ 
dictions and the measurements on the synthetic catalogues in the quasi-linear regime. 

Key words: cosmology: theory - large-scale structure of Universe - methods: ana¬ 
lytical - methods: statistical 


1 INTRODUCTION 

Galaxy clustering observations are a wealthy source of cos¬ 
mological information. The three-dimensional distribution of 
galaxies on large scales is usually characterized by means of 
two-point statistics such as the power spectrum (PS) or the 
two-point correlation function (2PCF). Galaxy redshift sur¬ 
veys such as the 2dF Galaxy Redshift Survey (Colless et al. 
2001), Sloan Digital Sky Survey (SDSS) I and II (York et al. 
2000), and the Baryonic Oscillation Spectroscopic Survey 
(BOSS, part of SDSS-III; Dawson et al. 2013) have demon¬ 
strated the significance that these probes have for preci¬ 
sion cosmology. Spectroscopic surveys that are ongoing or 
starting soon are the extended Baryonic Oscillation Spec¬ 
troscopic Survey (eBOSS, part of SDSS-IV; Dawson et al. 
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2015) and the Hobby Eberly Telescope Dark Energy Ex¬ 
periment (Hill et al. 2008). In the future, next-generation 
experiments such as the Dark Energy Spectroscopic Instru¬ 
ment (Levi et al. 2013), 4MOST (de Jong et al. 2014), the 
Subaru Prime Focus Spectrograph (Ellis et al. 2014) and the 
space-based Euclid mission (Laureijs et al. 2011) will allow 
for an even more accurate exploration of galaxy clustering 
over a wide range of redshifts, giving access to invaluable 
cosmological information. 

The strength of the clustering signal along the line-of- 
sight direction (LOS) differs from that in the transverse di¬ 
rections because of the peculiar velocities of the galaxies. 
This effect, dubbed redshift-space distortions (RSD; Kaiser 
1987), is mostly due to the large-scale infall of galaxies into 
gravitational potential wells and the non-linear velocity dis¬ 
persion due to the orbital motion of galaxies in dense, col¬ 
lapsed structures. By assuming a fiducial cosmology to con¬ 
vert redshifts to distances, measurements of the anisotropic 
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clustering pattern caused by RSD provide information on 
the linear growth of structure (Sanchez et al. 2013; Beutler 
et al. 2014; Samushia et al. 2014; Reid et al. 2014). Fur¬ 
ther, differences between the true cosmology of the Universe 
and the assumed fiducial one affect anisotropic clustering 
measurements through the so-called Alcock-Paczynski effect 
(Alcock & Paczynski 1979), which can be used in combina¬ 
tion with the signal from baryonic acoustic oscillations as 
a probe of the expansion history of the universe. Modelling 
this effect allows for measurements of the angular diameter 
distance Da{z) to and the Hubble parameter H(z) at the 
mean redshift z of the survey (see e.g. Anderson et al. 2014). 

A crucial ingredient of these cosmological analyses is 
an accurate model of the observed clustering statistics at 
cosmological, i.e., quasi-linear, scales, including the effect of 
non-linear evolution, galaxy bias and RSD. The validation 
of such models requires precise estimates of the covariance 
matrix of the measured anisotropic power spectrum or two- 
point correlation function. While the estimate of the cluster¬ 
ing covariance matrix for the galaxy sample itself can be gen¬ 
erated by brute-force production of hundreds or thousands 
of synthetic realizations using fast approximate schemes for 
the formation of cosmological structures, model testing is 
better performed with the full non-linear covariance of the 
clustering statistic. Hence, the model verification usually re¬ 
lies on large-volume N-body simulations of which only a 
few realizations can be produced due to limited run-time 
and memory. Covariance estimates from the ‘mock’ method 
suffer from the noise due to the finite number of realiza¬ 
tions affecting the inverse of the covariance matrix, which 
propagates into the obtained constraints (Taylor et al. 2013; 
Dodelson & Schneider 2013; Percival et al. 2014; Taylor & 
Joachimi 2014). The requirement to keep these error contri¬ 
butions small means that the number of mock realizations 
could reach the regime of a couple of thousands. The com¬ 
plexity of the problem increases further if the dependence of 
the covariance on the cosmological model is to be analysed 
as well (e.g., for the BAO covariance, Labatie et al. 2012). 
Then, the number of mock realisations significantly increases 
and the dependence of the covariance on any cosmological 
parameter can only be interpolated. 

Alternative methods can be applied to reduce the re¬ 
quired number of independent realizations. Hamilton et al. 
(2006) uses a reshuffling scheme for the phases of a peri¬ 
odic simulation to estimate the PS covariance from a sin¬ 
gle realization. A similar resampling was used by Schneider 
et al. (2011). Parametric models of the covariance matrix 
that only depend on a few free parameters would also be 
calibrated against a smaller number of simulations (Pear¬ 
son & Samushia 2015). The noise in the obtained covariance 
matrix can also be reduced by means of a shrinkage estima¬ 
tion (Pope & Szapudi 2008) or covariance tapering (Paz & 
Sanchez 2015). 

In this work, we present the most simplistic theoretical 
framework for the modelling of anisotropic clustering covari¬ 
ance. Previous studies showed that a Gaussian likelihood 
for the two-point statistics is sufficient for RSD experiments 
in the linear and quasi-linear regime (Manera et al. 2012). 
This paper aims at an extension of such covariance mod¬ 
elling for anisotropic clustering statistics in order to allow 
RSD analysis and model verification using synthetic cata¬ 
logues with periodic boundary conditions. So far, the the¬ 


oretical modelling of the covariance of anisotropic cluster¬ 
ing statistics has not been studied in detail - in contrast 
to the covariance of the monopole of 3D clustering statis¬ 
tics (Feldman et al. 1994; Smith et al. 2008; Sanchez et al. 
2008) or the covariance of angular 2-point statistics (Crocce 
et al. 2011). In this paper we present the explicit formu¬ 
lae for the covariance of the multipoles as well as clustering 
wedges of the redshift-space galaxy PS and 2PCF for cu¬ 
bic mock catalogues. For RSD and BAO studies, usually 
distance scales between 30 h~^ Mpc and 180 h~^ Mpc and 
wavenumbers between 0.02 hMpc“^ and 0.2 /iMpc“^ are 
considered as relevant (quasi-linear regime of gravitational 
evolution, BAO features within this range). Our model is 
verihed on these scales using a set of large-volume N-body 
simulation whose haloes have been populated with galaxies 
mimicking the clustering statistics of real surveys by use of 
the halo occupation distribution (HOD; Peacock & Smith 
2000; Scoccimarro et al. 2001; Berlind & Weinberg 2002) 
technique. The incorporation of the effects of a non-trivial 
survey geometry (described by the angular and radial selec¬ 
tion of galaxies) is left for future work. 

This paper is organized as follows: the anisotropic 
galaxy clustering statistics are introduced in Section 2.1. 
Measurements of their anisotropic moments and their co- 
variance are presented in Section 2.2. Section 2.3 presents 
our theoretical modelling of the anisotropic clustering covari¬ 
ance and gives the main results of this paper. In Section 3 
we validate our modelling using a set of large-scale N-body 
simulations which are described in detail in section 3.1. We 
show that our model successfully gives a smooth represen¬ 
tation of the data covariance of the simulations anisotropic 
clustering in section 3.3. Section 4 concludes our work. 


2 THE METHODOLOGY 

This section contains a description of the galaxy two-point 
statistics and their measurements. Considering the covari¬ 
ance of the two-point statistics, we describe how they are 
estimated from data and give a theoretical model for the 
Gaussian covariance matrix which is the main result of this 
paper. 

2.1 Galaxy two-point statistics 

Galaxy clustering measurements are based on the galaxy 
overdensity field. 


where n(x) denotes the actual number of galaxies around a 
point X, and n{x) is the expected galaxy number density. 
Most clustering analyses focus on the power spectrum (PS) 
and its Fourier transform, the two-point correlation function 
(2PCF). The 3D power spectrum P{k) gives the covariance 
of the overdensity field in Fourier space, 

{S{k) r{k')) = (27r)®fo(fc - k') [P{k) + h-i] , (2) 

where Sn{k — k') is the 3D Dirac delta function and S*{k) de¬ 
notes the complex conjugate of the Fourier transform of the 
overdensity, 5{k). For this relation we took into account that 
the discrete galaxy positions of a given survey are sampled 
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by a Poisson point process leading to a shot-noise contribu¬ 
tion n~^ (Feldman et al. 1994). Redshift-space distortions 
depend on the line-of-sight (LOS) parameter jj, which is de¬ 
fined as the cosine of the angle between the LOS and the 
separation vector of a galaxy pair. None the less, P{k) is 
still statistically invariant under rotations around the LOS. 
Hence, we express the power spectrum in terms of the abso¬ 
lute wavenumber and the LOS parameter, P{k) — P{k,fi). 

In order to have measurements with a significant signal- 
to-noise ratio, the anisotropic PS is usually projected onto 
multipole moments (Padmanabhan & White 2008), 

2 / 4-1 

P,{k) = J ^ Ctiu) P{k, u) d/r, (3) 


where Ci{u) denotes the Legendre polynomial of order 1. 
On the other hand, clustering wedges, which have first been 
defined in configuration space by Kazin et al. (2012), are 
projections into wide bins of the LOS parameter, 

Plt^ (fc) = - / P(fc, p) d/r. (4) 

- Ml Jmi 

where /ii and /12 define the lower and upper limit, respec¬ 
tively, of non-intersecting wedges such that /i2 — Mi — 

The total number of wedges defines Afi. We use the usual 
convention that in the case of only two wedges, these are 
labelled P± and Py for the M-ranges [0,0.5] and [0.5,1], re¬ 
spectively. 

Anisotropic clustering is often analysed in terms of the 
2PCF defined as ^(s) = {S{x) S{x + s)). Hence, it is the 
configuration-space counterpart of the PS, 

The 2PCF can also be decomposed in terms of the Legen¬ 
dre multipoles (Hamilton 1997); the decomposition can be 
derived from the PS multipoles using 


Cr(s) 




Pi(k)ji{ks) 


k^ dk, 


( 6 ) 


where je{x) is the spherical Bessel function of order £. Clus¬ 
tering wedges in configuration space are defined completely 
analogous to their Fourier space counterparts, 




1 


M2 — Ml 


J fll 


C(s,m) dM- 


( 7 ) 


By convention, we use the labels and for the case of 
two wedges. 

The relation between 2PCF wedges and multipoles is 


£ 


( 8 ) 


ber bins whose centres are at ki. The subscript x can ei¬ 
ther refer to a multipole index £ or to a clustering wedge 
limited by some {mi,M 2 }- Analogously, the 2PCF moments, 
have been measured in Nb separation bins cen¬ 
tred at Vi. 

Assuming that S{k) is a Gaussian random field, the 
^ 2 

Fourier modes |(5(fc)| follow a Rayleigh distribution (see e.g. 
Kalus et al. 2015). In the following, we assume that the num¬ 
ber of Fourier modes observed is large enough to validate the 
assumption that the power spectrum follows a multi-variate 
Gaussian distribution with fixed covariance.^ 

Given that the true inverse of covariance matrix, if), is 
known and fixed, we can compute the likelihood of a model 
prediction for the clustering statistic, P^ or , by^ 


^[pi I Pi,',p') 



1 

2 


xGP^P^ 



( 10 ) 


Here, the log-likelihood or function is given by 

x\n, {pi - Pi), ( 11 ) 

x,y i,j 


where 'ipxyij are the elements of the precision matrix il>. 

Usually the covariance matrix is estimated from a 
set of Njn mock measurements, denoted where n G 

{1,..., Nui}, so that® 

E - (P^)) - (Pi)) ■ (12) 

n 


Here, the mean over all simulations is given by 



n 


( 13 ) 


In order to account for the noise in (C^)“®, the likelihood 
of equation (10) must be marginalized over the distribution 
of these uncertainties which, in the case in which P^ are 
Gaussian random samples, follows an inverse Wishart dis¬ 
tribution with the true precision matrix r/;. In most analyses 
however, (C^)“® is treated as the exact inverse covariance, 
resulting in known biases of the parameter estimates which 
must be corrected for (see e.g. Kaufman 1967; Hartlap et al. 
2006; Percival et al. 2014; Taylor & Joachimi 2014). The 
smaller the number of realizations, Nm, the larger these cor¬ 
rections are; in the extreme case of Nh > Nui — 2, the co- 
variance matrix becomes singular. 

We want to avoid generating and processing a large 
number of synthetic catalogues and hence aim for a smooth 
and non-singular estimate of the covariance matrix by the¬ 
oretical modelling of the data covariance and its inverse. 


where is the average of a Legendre polynomial over the 
M-range of the wedge, 

I ^H4-A/j 

71^,4 = ^ y (9) 

This integral is evaluated in equation (A22) in the appendix. 

2.2 Measurements of the covariance matrix of 
two-point clustering measurements 

We consider that the moments of the anisotropic PS have 
been measured to be = Px{ki) in Nh different wavenum- 


^ Deviations from the Gaussianity of the angular power spectrum 
likelihood have recently been analysed (Sun et al. 2013; Kalus 
et al. 2015) and it was shown that constraints on primordial non- 
Gaussianity are affected by the negligence of the non-Gaussian 
contribution. The logarithm of the spatial power spectrum was 
found to be closer to a Gaussian random variable at large scales 
than the power spectrum itself by Ross et al. (2013). We neglect 
the dependence of the covariance on cosmological parameters as 
suggested by Carron (2013). 

^ For the 2PCF likelihood, replace the P^ and P^ by and 
^ The 2PCF covariance matrix is likewise estimated from 

mock measurements 
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2.3 Modelling of the covariance of galaxy 
two-point clustering measurements 

2.3.1 Fourier Space: the Power Spectrum Covariance 




2 {2-k)* 


■L 


ki+Ak/2 


-I^k/2 


Ak, (16) 


The assumption that the 2D power spectrum follows a Gaus¬ 
sian distribution leads to the following relation for the PS 
mode-by-mode covariance (see Feldman et al. 1994, for the 
monopole), 


Cov[P(fe),P(fc')] = ^-^^fo(fe-fc') [P(fc,/r)-fh-Y, 

Vs 

(14) 

where 14 is the volume of the sample."^ The Dirac delta func¬ 
tion reflects the independence of Fourier modes for a random 
held with statistical translational invariance and in the ab¬ 
sence of gravitational mode coupling. The anisotropy of the 
PS is taken into account by the LOS dependency of P(fc, p). 

Here, we neglected the trispectrum contribution (Scoc- 
cimarro et al. 1999) and the super-sample covariance (SSC; 
also called beat-coupling). Further, the fact that the local 
density estimate is obtained in the presence of power on 
scales larger than the survey size affects the power spec¬ 
trum covariance of real surveys in the same manner as SSC, 
but the net effect only makes up 10% of the original beat- 
coupling effect (de Putter et al. 2012). The contribution of 
the two latter effects to the PS covariance has been shown 
to become important at non-linear scales k > 0.1 /iMpc“^ 
(de Putter et al. 2012; Takada & Hu 2013; Li et al. 2014). 
The mode coupling in the standard trispectrum term due to 
non-linear gravitational evolution is subdominant compared 
to SSC, but adds to it creating a plateau in the signal-to- 
noise ratio of the matter power spectrum in the mildly non¬ 
linear regime and beyond (Carron et al. 2015). The mod¬ 
elling of these effects is beyond the scope of this paper and 
their analysis is left for future work. 

The theoretical covariance matrix for binned anisotropic 
PS measurements is obtained by averaging over the num¬ 
ber of independent Fourier modes k that contribute to each 
wavenumber bin. We assume the bins to be centred at ki 
with width Afc. For readability, the derivation of the follow¬ 
ing covariance expressions are presented in appendix A. A 
useful definition is the multipole expansion of the per-mode 
covariance given by 


2 N _ ( 2^1 + 1 ) ( 2^2 + 1 ) 


/: 


14 

1 r 


p(fc, p) — 
n 


Cej_{p) Ci^ip) dp. (15) 


The normalization was chosen such that the pre-factor of 
crQo(fc) is the usual 2/14 for an isotropic PS. Equation (A5) in 
appendix A provides the corresponding formula expressing 
the power spectrum in a multipole series. 

As shown in section Al, our ansatz for the bin-averaged 
PS multipole covariance {ki, kj) = Cov [PlwK] yields 


where the volume of the shell in fc-space is 14^ = 47r[(fci -f 
Afc/2)^ — {ki — Afc/2)^]/3. The full multipole expansion pre¬ 
sented in the appendix shows that the multipole covariance 
matrix has terms which mix all multipole contributions. This 
implies that the monopole covariance is not only given by 
terms depending on {Po{k) h“^)^, but also contains con¬ 
tributions from all higher-order multipoles as well. 

Second, we consider clustering wedges P/*/(A:) as de- 
Hned by equation (4). We denote the measurement of 
p^+An wavenumber bin centred around ki as P/. 

Then, the theoretical covariance of the bin-averaged wedges, 


C^^/(fci, kj) = Cov j^P/, P^,j, is given by (see section A2) 


,(k k ) — ^ J ■ S 

Kj 1 — y2 


^fcj+Afc/2 

1 ■ 

^ / / 

P{k,~p) + ^ 

J ki — A/c/2 J /i 

n 


dp k^ dk, 


(17) 


where <5^^/ reflects that the wedges are discrete and non¬ 
intersecting. In contrast to the multipoles, different Fourier 
space wedges are not correlated in our linear Gaussian the¬ 
ory, even after integration over the wavenumber bin. 


2.3.2 Configuration Space: the 2PCF Covariance 

As in the previous section, we only present here the final 
results and give the full derivations in section A3. By dehn- 
ing the per-mode multipole covariance as in equation (15), 
we can write the bin-averaged 2PCF multipole covariance 
matrix, (si, 8j) = Cov , as 

/ k^ o-ly.^{k)jti{ksi)ft.2{kSj)dk, 

(18) 

where je{ksi) is the bin-averaged spherical Bessel functions 
as defined in equation (A19). The bin average has been 
shown to be necessary in order not to overestimate the 2PCF 
covariance (Cohn 2006; Sanchez et al. 2008). 

In section A4, we derive the covariance of the 2PCF 
clustering wedges, Cf^y{si,Sj) = Cov to be 

^1,^2 

poo 

X / <^}ii2{k)jti{kSi)ft2{kSj)dk. (19) 
Jo 

Here, is the average of a Legendre polynomial over the 
p-range of the wedge as defined in equation (9). 


^ In our case of a periodic box, the volume is given by 14 = 
, where L is the side length of the box. For a survey with 
a selection function given by a random catalogue with varying 
number density, the effective volume has to be estimated from 
the window function (Bernstein 1994; de Putter et al. 2012). 


3 VALIDATION OF THE MODEL 

In this section, we first describe the set of N-body simula¬ 
tions that we will use to validate our model of the covariance 
of anisotropic galaxy clustering two-point statistics. Next, 
we show and discuss the agreement between the theoretical 
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Table 1. The cosmological parameters of our set of Minerva 
simulations. 


parameter ^ crs 

value 0.285 0.715 0.695 0.9632 0.828 


Table 2. The parameters of the HOD model of equations (20) 
to (22) defining our ‘CMASS-like’ galaxy sample at 2 : = 0.57. All 
masses are in units of h~^ M©* 


parameter 

value 

logio(Mmm) 

13.07 

logio(Mo) 

13.1 

logio(M{) 

14.2 

ClogiQ M 

0.347 

a 

0.8 


'e' 

-d 


10'“ 

10'’ H- 

lot I 

10; 

10"' 

10-; 

10 '^ E- 
10'® 

10’'“ 


SubFind haloes 
Tinker et al. 



prediction and the results from the simulations. We specif¬ 
ically address the question of whether the predictions pre¬ 
sented in the previous section are precise enough to allow 
performance test of RSD models for the clustering statistics 
of current galaxy surveys. 


3.1 The N-body Simulations and Galaxy 
Catalogues 

The Minerva simulations are a set of 100 N-body simula¬ 
tions run using Gadget® (last described in Springel 2005) 
with 1000® dark matter (DM) particles per realization in 
a cubic box of side length 1500 h~^ Mpc with periodic 
boundary conditions. The simulations were started at red- 
shift Zini = 63 using 2LPT initial conditions. The input lin¬ 
ear power spectrum for the initial conditions was calculated 
using Camb (Lewis et al. 2000) for the cosmological pa¬ 
rameters chosen for our set of simulations (see Table 1), 
which match the best-fitting ACDM model of the WMAP9 
+ BOSS DR9 ^(r) analysis (Sanchez et al. 2013, Table I). 
The positions and velocities of the DM particles were stored 
for five output redshifts « G {2.0,1.0,0.57,0.3,0}, jointly 
with the halo positions and velocities found by a friends-of- 
friend halo finder applied with the standard linking length 
equal to 0.2 of the mean inter-particle separation. In a post¬ 
processing step, SubFind (Springel et al. 2001, section 4.2) 
was run to generate the final halo catalogues for each realiza¬ 
tion and output redshift. The mean halo mass function for 
2 : = 0.57 is plotted in the upper part of Fig. 1. The minimum 
resolved halo mass is rrimin = 2.67- 10®® /i"®Mq. 

The volume of each realization is large enough to al¬ 
low for the analysis of anisotropic galaxy clustering probes 
with a precision comparable to present-day galaxy redshift 
surveys. In order to generate a galaxy sample comparable 
to the CMASS sample of BOSS, we populate the haloes 
and subhaloes of each simulation with galaxies according to 
a halo occupation distribution (HOD) model with suitable 
parameters. In such HOD models the average number N of 
synthetic galaxies in haloes of mass M is given by the mean 


® The latest public release is Gadget- 2 which is available at 
http://www.gadgetcode.org/. 


Figure 1. Upper panel: the mean halo mass function of the Sub¬ 
Find output for our Minerva snapshots at z = 0.57 (black points) 
and the prediction (gray dashed line) based on the recipe of Tin¬ 
ker et al. (2008). Lower panel: the mean galaxy occupation func¬ 
tion as defined in equation (20), {N(M)) (red solid line), and its 
decomposition into central (blue dashed line) and satellite com¬ 
ponents (green dotted line) for the parameters given in Table 2. 


occupation function N{M). We follow the parametrization 
of Zheng et al. (2007) and decompose this function into con¬ 
tributions from central and satellite galaxies, 

{N{M)) = (iVcen(M)) -F (W,t(M)), (20) 


where the mean central occupation function. 


(Vce„(M)) = 


1 + erf 


/ log M — log Mu 

\ G 'log A4 


( 21 ) 


has a smooth cut-off at Mmin modelled by an error function 
with relative scale uiogM to describe the scatter between 
the galaxy luminosity and mass. By setting the satellites 
contribution to 

(Wat(M)) = (Wen(M)) ( ^ ^ ^ , (22) 

we only assign satellites to haloes that are already populated 
by a central galaxy. The satellite galaxies are sampled with 
a Poisson distribution with the mean set by equation (22) 
where the cut-off mass scale is Mo, the normalization mass 
scale is M[, and the power-law slope is a. If a central galaxy 
is assigned to a halo, its position and velocity are derived 
from the most-bounded DM particle in that halo. The po¬ 
sitions and velocities of satellite galaxies are drawn from 
random DM particles associated to the hosting halo. 

The parameters chosen for our HOD sample at redshift 
2 = 0.57, the mean redshift of the BOSS CMASS sample, 
are given in Table 2; the resulting mean galaxy occupation 
function is shown in the lower panel of Fig. 1. The final 
synthetic galaxy catalogues have a mean galaxy density of 
h « 4 • 10“'* /i®Mpc“®. By comparing the redshift-space 
2PCF monopole of the HOD sample with the real-space one 
of the DM particles, we find a linear bias of 5® = 4.02 de¬ 
rived from the ratio at pair separations 40 /i“* Mpc < s < 
60 /i“® Mpc. Using the ratio of the real- and redshift-space 
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s [h-^ Mpc] 


Figure 2. Comparison of the 2PCF clustering wedges 
and of our HOD sample (red and blue solid lines, re- 

spectively, standard deviation over our 100 realizations indicated 
by the filled region) compared with the corresponding measure¬ 
ments from the CMASS sample of BOSS DRll (transverse wedge: 
circles, parallel wedge: diamonds) by Sanchez et al. (2014). 

monopole of the HOD sample at these separations, we find 
that the Kaiser factor is S' = 1.28, in perfect agreement with 
the theoretical value ofS = 1 + + |/3^ (Kaiser 1987) 

derived from the growth factor / = 0.76 for the Minerva 
cosmology. Here, we used p = f /b, where / is the growth 
rate, i.e., the logarithmic derivative of the growth function 
w.r.t. the scale factor. 

The agreement between the anisotropic clustering of the 
HOD galaxy sample and the CMASS sample of BOSS DRll 
can be seen in Fig. 2, where we show the comparison of the 
2PCF clustering wedges of both samples. 

No super-survey modes have been taken into account 
for our simulations (e.g. by use of the separate universe re¬ 
sponse, Li et al. 2014) in agreement with the absence of 
those modes in our modelling of the PS covariance in equa¬ 
tion (14). 

3.2 The two-dimensional power spectrum 

As the simplest case, we assume a linear prediction for the 
2D galaxy power spectrum where redshift-space anisotropies 
are caused by the linear Kaiser effect (Kaiser 1987), 

P“’'(A:,m) = &"PL(A:)(l-b/3V")", (23) 

where Ph{k) is the linear matter power spectrum, and b is 
the linear galaxy bias. 

With these assumptions, the only non-vanishing PS 
multipoles are the monopole, quadrupole, and hexadecapole. 

Fig 3 shows the mean PS monopole, quadrupole and 
hexadecapole from the Minerva simulations. The linear the¬ 
ory definition of equation (23), shown by the dotted lines, 
gives an inadequate description of the anisotropic galaxy 
power spectrum in the quasi-linear regime. To improve upon 
this description we perform a smoothing spline interpolation 
of our mean PS multipole measurements, Pf™°°*’'(fc), shown 



Figure 3. Comparison of the mean power spectrum multipoles 
of our HOD sample (points) with the linear-theory predictions, 
Pj*'“(fc), of equation (23) (dotted lines) and the smoothed interpo¬ 
lation, P™°°**'(fc) (dashed lines). Both the linear and smoothed 
multipoles are taken as input for linear and non-linear predictions 
of the Gaussian covariance. 

by the dashed lines in Fig 3, to create a noiseless non-linear 
power spectrum that we can use for covariance predictions, 

P“‘(fe,M) = Po™°°“'(A:)-bP|”°°“'(A:)/:2(M)+Pr°°*‘'(fc)T4(At) 

(24) 

For the estimation of the smoothing length, we take the 
measured dispersion of the PS multipoles into account. The 
BAO wiggles in the smoothed quadrupole have been slightly 
damped by this procedure but, due to the small signal-to- 
noise ratio of the BAO feature in the quadrupole, this does 
not affect the predicted covariance. 

In the same way as the linear model of equation (23), 
by definition this ansatz also has only monopole, quadrupole 
and hexadecapole contributions. Higher-order multipoles are 
negligible for k < 0.25 hMpc“^. In the following, we will use 
these two power spectrum models for the analysis and val¬ 
idation of our Gaussian covariance predictions (referred to 
as ‘lin’ and ‘smooth’, respectively). In the case of an appli¬ 
cation of the Gaussian covariance model for the use of RSD 
model performance tests, the smoothed power spectrum can 
be replaced by a preliminary fit of the RSD model to the 
data. 

3.3 Validation of the Accuracy of the Gaussian 
Covariance Model 

In order to compare the theoretical predictions presented in 
section 2.3 with the noisy estimates from the Minerva sim¬ 
ulations, we chose two different binning configurations: one 
leading to an invertible covariance matrix (in case of mul¬ 
tipoles up to the hexadecapole or three clustering wedges) 
denoted ‘large’ and one appropriate for fitting of CMASS- 
like measurements (for which the data covariance matrix 
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Table 3. The binning choices for the validation of our covariance 
model. 


space 

range 

A small 
A large 

unit 

number of bins 

s 

0 - 180 

5 

h~^ Mpc 

36 



15 


12 

k 

0 - 0.25 

0.005 

hMpc-i 

50 



0.010 


25 


obtained from 100 catalogues is singular) denoted ‘small’. 
These setups are listed in Table 3. 

The measurements of the redshift-space two-point clus¬ 
tering assume a plane-parallel LOS along one of the axis 
of the simulation box (‘distant-observer approximation’). In 
order to reduce the level of noise in the mean and covariance 
of the measurements, we measure the two-point statistic of 
every realization by assuming the LOS to be parallel to each 
of the three different axes and then average the results. 

Due to this averaging over the three LOS axes, the error 
on the covariance cannot easily be predicted analytically. 
Thus, we measure the error on the covariance matrix as 
jackknife estimate from the same set of mock measurements, 

(ACf,,,)' = ^ ' , (25) 

m 

where is the covariance estimate based on leaving 

out the mth realization, 

n^m 

(26) 

The error on the covariance matrix of all configuration space 
measurements, are obtained in an analogous way. We 

checked that the jackknife estimate of the error on the co- 
variance matrix for the case of a single LOS axis is in close 
agreement with the theoretical prediction given in Taylor 
et al. (2013). 

The calculation of the likelihood in equation (10) re¬ 
quires the estimation of the inverse of the covariance ma¬ 
trix, the precision matrix. If this matrix is taken to be the 
inverse of the data covariance, it is biased due to the sample 
variance. The level of noise is quite high in our case due to 
the small number of realizations and large number of mea¬ 
surements bins. Since for many analysis configurations, the 
data covariance is even singular, we do not comment on the 
accuracy of the precision matrix here. We shortly discuss the 
inverse data covariance matrix from our set of simulations 
for one particular binning scheme in appendix B. 

3.3.1 Fourier Space Covariance 

The predicted covariance of the PS multipoles is compared 
against the data covariance measured from our 100 simula¬ 
tions in Fig. 4 for the small (upper panel) and large (lower 
panel) binning setup. We show the dispersion of PS mul¬ 
tipoles given by {ki) = [Cu{ki, . The error bar 

is likewise given by Aaf{ki) = [ACff{ki,ki)Y^'^, where 
ACf({ki, ki) is estimated with the jackknife estimator given 




Figure 4. Comparison of the dispersion of the PS multipoles, 
of our HOD realizations (points) with the Gaussian pre¬ 
dictions from equation (16) for the small (upper panel) and large 
(lower panel) binning schemes. The dotted lines were derived us¬ 
ing the linear model of equation (23), while the dashed lines cor¬ 
respond to the smoothed non-linear recipe of equation (24). 

in equation (25). The difference between linear and non¬ 
linear predictions is small for the monopole and quadrupole 
dispersion. However, the linear prediction of the hexade- 
capole dispersion slightly overpredicts (underpredicts) the 
actual dispersion in the data for lower (higher) k, respec¬ 
tively, which is better matched by the results obtained from 
the smoothed multipole measurements. In appendix Al, we 
comment on the increase of the higher-order multipole vari¬ 
ance showing predictions ap^ and apg. 

The measurements of the off-diagonal terms 
Cf^ f^{ki,kj) (for £i ^ (.2 and/or ki yf kj) suffer from 
low signal-to-noise and hence show larger relative scat¬ 
ter than the dispersion terms. As the differences in our 
linear and non-linear modelling are much smaller than 
the variations in the measurements we find equally good 
agreement between the two predictions and our simulation 
data. Because of the scatter, we cannot identify additional 
contributions such as non-gaussianities in our covariance 
measurements. To highlight the full covariance properties, 
we define the correlation matrix for PS multipoles as 

Re.^t2(G,kj) = Cij^i^{ki,kj) l^cTpni(fei) CTpju (fej)j . (27) 

The upper panel of Fig. 5 we show a color representation of 
the data correlation matrix (upper triangular part) next to 
the theoretical prediction (lower triangular part). The lower 
panel shows a cut through the matrix. Note that we nor¬ 
malized the data covariance by the theoretical dispersion 
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fcj \h Mpc ^ ] 



Figure 5. Upper panel: The full correlation matrix of the PS 
multipoles, = Cf^^^(ki,kj)[a p^\{ki) a 

(normalized by the theoretical prediction), shows a dominant 
diagonal and a significant cross-correlation between monopole 
and quadrupole (shown here: large binning). Lower panel: Cut 
through the correlation matrix for = 0 at = 0.095 /iMpc“^. 
The correlation contamination from physical effects not accounted 
for by our model and from noise is at the 15% level and well within 
the error bars. 


obtained from the smooth non-linear power spectrum, Opw , 
so that there is no further noise contamination. The level 
of unaccounted data correlation in the off-diagonal terms is 
up to 20%, but this is most likely due to noise contami¬ 
nation instead of systematics not included in our Gaussian 
ansatz. The noise level for the monopole is larger than for 
the higher-order multipoles because of the average over the 
LOS directions which reduces the noise of the RSD but not 
of the real-space matter clustering. 

Assuming three PS wedges, we label the covariance be¬ 
tween the n-th and m-th wedge as kj) and the 

dispersion as for clarity. For 

the case of two clustering wedges, we stick with the labels 
T and ||. Fig. 6 compares the predicted covariance of two 
(upper panel) and three (lower panel) wedges to the data 
covariance measured from our 100 realizations for the small 
binning configuration. The same comparison is plotted for 
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Figure 6 . Comparison of the dispersion of the PS wedges of 
our HOD realizations (points) for the small binning configuration 
with the Gaussian predictions from equation (17) for the linear 
(dotted lines) and smoothed non-linear (dashed lines) input power 
spectra. The upper and lower panels show the cases of two wedges, 
o-p^ II three wedges, o-p^^ respectively. 



[h Mpc ^ ] 


Figure 7. Comparison of the rms of the PS wedges of our HOD 
realizations (points) for the large binning configuration with the 
Gaussian predictions from equation (17) for the linear (dotted 
lines) and smoothed non-linear (dashed lines) input power spec¬ 
tra. The upper and lower panels show the cases of two wedges, 
ap^ II three wedges, o'p^^ respectively. 
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fcj [h Mpc ^ ] 





Sj [h ^ Mpc] 


Figure 9. Comparison of the dispersion of the 2PCF multipoles, 
(si) of our HOD realizations (points) compared with the Gaus¬ 
sian predictions from equation (18) for the small (upper panel) 
and large (lower panel) binning configuration. The dotted and 
dashed lines correspond to the linear and smoothed non-linear 
input power spectra, respectively. In addition, we show the predic¬ 
tion of the real-space covariance as given in Sanchez et al. (2008), 
rescaled by the linear Kaiser factor S = 1.28 as dot-dashed black 
line. This curve cannot be distinguished from the Gaussian pre¬ 
diction for the monopole given in this work. 


Figure 8. Upper panel: The full correlation 

matrix of the PS wedges, -R3w,n,m(^i? ^j) = 

(fc*)o'pni , shows a domi- 

nant diagonal and no significant cross-correlation between 
wedges (shown here: large binning scheme). Lower panel: Cut 
through the correlation matrix for the most transverse wedge, 
?3w,li a-t ki = 0.095 /iMpc“^. The correlation contamination 
from physical effects not accounted for by our model and from 
noise is at the 15% level (same level as the error bars). 

the case of the large bins in Fig. 7. In the linear-theory pre¬ 
dictions, the BAO peaks and troughs in the dispersion are 
not as damped as in the data, especially for the most parallel 
wedges. Further, using linear theory the high-fc dispersion of 
the intermediate and most perpendicular wedge is underpre¬ 
dicted (fc > 0.15 h/Mpc). 

We complete the PS wedge covariance analysis by show¬ 
ing the full correlation matrix compared to the non-linear 
prediction in Fig. 8. Here we restrict the analysis to the 
3-wedges case for brevity; the 2-wedges correlation ma¬ 
trix has similar properties. As described in section 2.3, 
the cross-correlation fcj) vanishes for n ^ m or 

ki ^ kj. The data correlation matrix, R3vr,n,m{ki, kj) = 
Cil,^n,m{ki,kj)[apn^ ^{ki)apn^ ^{kj)]~^, again normalized 
by the non-linear theoretical prediction, shows a level of un¬ 
accounted cross-correlation of up to 15%. As for the multi¬ 
poles, this is most likely noise contamination and the scatter 


cannot be associated with any systematic deviations from 
our Gaussian model. 

3.3.2 Configuration Space Covariance 

In configuration space, we define the dispersion of 2PCF 
multipoles by cr|(si) = [C|^(si, Si)]^/^. In Fig. 9, the pre¬ 
dicted dispersion of the 2PCF multipoles is compared to the 
data dispersion measured from our 100 simulations for the 
small (upper panel) and large (lower panel) binning setup. 
The difference between linear and non-linear predictions for 
the multipole covariance is negligible compared to the devi¬ 
ations between the measurements and the theory estimates. 
Those deviations are largely due to the low number of real¬ 
izations and the fact that the Fourier transform to configu¬ 
ration space translates uncorrelated modes into highly cor¬ 
related measurements of the binned anisotropic 2PCF. The 
dispersion of higher-order multipole moments of the corre¬ 
lation function is discussed in appendix A3. 

The correlation between measurement bins is best il¬ 
lustrated by the full information encoded in the correlation 
matrix, which is defined in an analogous way to the case of 
PS multipoles, 

-R|ir2 ~ *^5“^ (®j)j ■ (28) 

Cross-covariance terms, (si, Sj), show a large scatter 
due to the low signal-to-noise of our measurement. We find 
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Figure 10. Upper panel: In the full correlation 
matrix of the 2PCF multipoles, — 

Cp p (si, Sj)[a^ni (si) {sj)]~^ (normalized by the theoretical 

prediction), shown here for the large binning, we see a high level 
of cross-correlation which is only slowly declining away from the 
main diagonals and the diagonals of the monopole-quadrupole 
and quadrupole-hexadecapole sub-matrices. The monopole and 
the hexadecapole are not correlated. Lower panel: Cut through 
the correlation matrix for £i = 0 at = 112.5 h~^ Mpc. 
Although the correlation curves inferred from our simulations 
are noisy, they follow the same trends as out our theoretical 
predictions. 


that the differences in our linear and non-linear modelling 
are much smaller than the scatter in the measurements and 
we conclude agreement between the different predictions and 
data on an equally good level. 

In the upper panel of Fig. 10, we show the structure 
in the data correlation matrix (upper triangular part) next 
to the theoretical prediction (lower triangular part). For 
a better visualization of the level of correlation, the lower 
panel shows a cut through the matrix. The prediction for 
the monopole, quadrupole and hexadecapole correlation is 
in very good agreement with the measurements. We see how¬ 
ever that the data cross-correlations between different mul¬ 
tipole moments (especially between monopole and hexade¬ 
capole) are noisy. 



Figure 11. Comparison of the dispersion of the 2PCF wedges of 
our HOD realizations (points) for the small binning configuration 
with the Gaussian predictions from equation (19). The dotted and 
dashed lines correspond to the linear and smoothed non-linear in¬ 
put power spectra, respectively. The upper and lower panels show 
the cases of two wedges, n (si), and three wedges, 
respectively. 



s. [ft-iMpc] 


Figure 12. Comparison of the dispersion of the 2PCF wedges of 
our HOD realizations (points) for the large binning configuration 
with the Gaussian predictions from equation (19). The dotted and 
dashed lines correspond to the linear and smoothed non-linear in¬ 
put power spectra, respectively. The upper and lower panels show 
the cases of two wedges, || (si), and three wedges, 
respectively. 
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As for PS wedges, we label the covariance between the 
n-th and m-th 2PCF wedge as S'lid the dis¬ 
persion as = [Clw.n, 71 ( 5*1 In the case of two 

wedges only, we label each one with a subscript _L or |{. 

In Fig. 11, we plot the predicted dispersion of two (up¬ 
per panel) and three (lower panel) wedge measurements in 
configuration space compared to the data dispersion mea¬ 
sured from the Minerva simulation for the small binning 
configuration. These predictions are based on equation (19), 
including contributions up to £ 1 , £2 < 6. We discuss the con¬ 
vergence of the wedge dispersion with the inclusion of higher- 
order multipoles in appendix A4. The same comparison is 
plotted for the case of the large bins in Fig. 12. In both cases, 
there is no significant difference between the linear and the 
non-linear predictions. The estimate from linear theory is 
always slightly larger than the estimate from the smoother 
interpolation of the measured power spectrum multipoles, 
but this difference is much smaller than the deviation be¬ 
tween theory and data. 

The cross-correlation for 2PCF wedges is as signifi¬ 
cant as for the configuration-space multipoles. In analogy 
to the discussion on the Fourier space wedges, we define 
the correlation matrix for 2PCF wedges as n(5*,5j) = 
CL,n,m(5*.5l)kej4,_^(50cr5n4,^(sj)]“\ normalized by the 
non-linear theoretical prediction. We show the full correla¬ 
tion matrix compared to the non-linear prediction in Fig. 13. 
Again, we restrict our discussion to the case of three wedges 
for brevity. The correlation matrix for two 2PCF wedges has 
similar properties. As the Fourier transformation mixes inde¬ 
pendent Fourier modes, the covariance C'|w,n,m(5»! 
decays slowly with increased separation between the dis¬ 
tance bins and ranges of the LOS parameter. The data corre¬ 
lation (upper panel of Fig. 13) shows good overall agreement 
between the theoretical prediction (lower triangular part) 
and our measurement (upper triangular part). The plotted 
full matrix shows also that the cross-correlation structure 
for 2PCF wedges is more complex than for 2PCF multipoles. 
For better visualization, we show a cut through the matrix 
in the lower panel of Fig. 13. We find differences between 
our measurements and the theoretical predictions of up to 
20%, which is a bit higher than for the multipole case. As for 
the previous cases, these discrepancies are most likely due 
to noise contamination and the level of noise is very similar 
for different cuts through the matrix. The model under- or 
overpredicts the data over wide ranges of pair separations 
because the noise in the covariance is as correlated as the 
2PCF itself. Hence, we are not able to test for deviations 
from of our Gaussian ansatz with the number of realizations 
at hand. 

As already discussed in section 2.3, the shortcomings of 
our model for the analysed sample of N-body simulations 
is the neglected trispectrum contribution which is only rel¬ 
evant for scales in the non-linear regime. We would need a 
much larger set of N-body simulations in order to validate 
extensions to the Gaussian model. 

As a final test, we analyse whether our model for the 
monopole covariance has improved over the simplifying as¬ 
sumptions used by Sanchez et al. (2008), where the Gaussian 
real-space covariance from linear theory has been rescaled 
by the Kaiser factor for linear RSD. This model defines the 




[h ^ Mpc] 


Figure 13. Upper panel: The full correlation matrix of the 
2PCF wedges, Corr[j 3 „^i(fci),^ 3 w_„(fcj)], is dominated by a high 
level of correlation between different wedges and distance bins 
(shown here for the large binning configuration). Lower panel: 
Cut through the correlation matrix for the most transverse 
wedge, ^ 3 w.l, at s; = 112.5 /i“^Mpc. Although the data cross¬ 
correlations are noisy, they are well described by our theoreti¬ 
cal prediction (the significant correlation of the covariance is dis¬ 
cussed in the text). 


monopole covariance in real-space as follows: 


^S.lin.rs/ 

'-* 0,0 15 *> 



b^PGk) + - 

n 

X Jt{ksi)ji{ksj) k^ dk, 


(29) 


From this, the a-space covariance is estimated to be 
s,) = Sj), where 5 = 1 + §/?+ 

As shown in Fig. 9, this simple ansatz gives already an excel¬ 
lent description of the covariance of the redshift-space 2PCF 
monopole, even though the anisotropy of the input power 
spectrum has been neglected. Any RSD or BAG fit using 
anisotropic clustering measurements must take the higher- 
order multipoles into account which can only be done by the 
formulae presented in this work. 

Additional tests indicate that the discrepancy between 
the simplified monopole recipe and the full Gaussian model 
presented here depends on the bias. While there is a neg- 
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ligible difference for the highly biased HOD galaxy sample 
analysed here (where = 4.02), the covariance of an nnder- 
biased tracer sample (& J; 1, similar h) is only correctly mod¬ 
elled taking the full anisotropy of the input power spectrum 
into acconnt. 


4 CONCLUSIONS 

In this work, we presented explicit formulae for the Gaussian 
covariance matrix of anisotropic two-point clustering mea¬ 
surements. Aiming at precise covariance estimates for the 
verification of models of the redshift-space clustering two- 
point statistics, we looked at the covariance of clustering 
wedges (i.e., large bins in the line-of-sight parameter) and 
multipole moments of the clustering signal, both in con¬ 
figuration space as well as in Fourier space. The formulae 
presented here rely on a model for the input power spec¬ 
trum. In our analysis, we tested two different prescriptions 
for P{k, p): first, a linear model for the redshift-space galaxy 
power spectrum based on the linear Kaiser effect and a lin¬ 
ear galaxy bias on top of the linear-theory prediction for the 
matter power spectrum. Second, we use a smoothed inter¬ 
polation of the measured galaxy power spectrum from our 
simulations. 

The covariance model has been tested and validated 
against a set of large-volume N-body simulations. The num¬ 
ber of realizations is enough to have a precise mean mea¬ 
surement of the non-linear anisotropic clustering signal, but 
the recovered covariance matrices exhibit significant noise 
due to the relatively small number of realizations. 

We show that the Gaussian covariance from the inter¬ 
polated non-linear power spectrum describes the measured 
clustering covariance very well. Not only is the overall shape 
of the dispersion of PS wedges and multipoles in excellent 
agreement, also broad features such as damped BAOs are 
described accurately. Relying on the linear power spectrum 
model, the Gaussian prediction is also in very good agree¬ 
ment with the observed data covariance, but the BAG fea¬ 
tures in the PS covariance are slightly overpredicted. In con¬ 
figuration space, the measurements of the anisotropic clus¬ 
tering are highly correlated. The theory predictions for the 
covariance of wedges and multipoles of the two-point correla¬ 
tion function (2PCF) accurately describe the full correlation 
structure. The differences in these covariances resulting from 
the two different descriptions for the input power spectrum 
are almost negligible. 

The aim of this work is limited to predictions of the 
anisotropic clustering covariance on quasi-linear scales in 
order to allow RSD and BAG fits of the clustering mea¬ 
surements from a limited set of large-volume N-body real¬ 
izations. Further work is needed to incorporate the beat¬ 
coupling with super-survey modes and the contributions 
from the connected trispectrum part of the covariance with 
arises due to non-linear evolution (Scoccimarro et al. 1999). 
These effects have been neglected here due to absence 
of super-survey modes (super sample covariance; SSC) in 
the N-body realizations and the restrictions to quasi-linear 
scales {k < 0.2 h,Mpc“^, s > 40 h~^ Mpc). In order to vali¬ 
date a model for the influence of SSG on the anisotropic clus¬ 
tering covariance, an similar analysis to the work in de Put¬ 
ter et al. (2012); Takada & Hu (2013); Li et al. (2014) - 


i.e., taking account of the local density estimate and the 
beat-coupling of smaller scales with super-survey modes - is 
required, which is beyond the scope of the work presented 
here. 

We expect careful modelling to be more important the 
more complex the analysed galaxy sample is, and the deeper 
in the non-linear regime the clustering probe advances. The 
Gaussian model derived here can easily be extended to ac¬ 
count for the dependence of the covariance on cosmological 
parameters by varying the input power spectrum. In future 
work, we plan to extend our studies to models of the cluster¬ 
ing covariance of surveys with non-trivial geometry - for the 
analytical treatment in case of isotropic clustering measure¬ 
ments, see de Putter et al. (2012). After submission of our 
initial manuscript, G’Connell et al. (2015) published a pre¬ 
print presenting a way to take into account non-Gaussian 
contributions to the anisotropic 2PGF covariance matrix for 
a CMASS-like survey. 

As final remark, we want to emphasize that the set 
of Minerva simulations is especially suited for RSD model 
testing in the mildly non-linear regime and such tests can 
make use of the covariance matrices presented in this work. 
Among other tests, we will apply this in the context of the 
RSD analysis of the final release of the BGSS galaxy clus¬ 
tering sample currently being prepared. 


ACKNOWLEDGEMENTS 

We thank Daniel Farrow, Elisabeth Krause, Martha Lip- 
pich, and Francesco Montesano for helpful discussions. We 
would like to thank the anonymous reviewer of our ini¬ 
tial submission for many helpful comments. JNG, AGS and 
SSA acknowledge support from the Transregional Gollab- 
orative Research Centre TR33 ‘The Dark Universe’ of the 
German Research Foundation (DFG). CDV acknowledges 
support from the Ministry of Economy and Competitive¬ 
ness (MINECG, Spain) through grants AYA2013-46886 and 
AYA2014-58308. The analysis has been performed on MPE’s 
computing cluster for the Euclid project and the Hydra clus¬ 
ter at the Rechenzentrum Garching (RZG). 


References 

Alcock C., Paczynski B., 1979, Nature, 281, 358 
Anderson L., et al., 2014, Mon.Not.Roy.Astron.Soc., 441, 24 
Berlind A. A., Weinberg D. H., 2002, Astrophys.J., 575, 587 
Bernstein G. M., 1994, Astrophys. J., 424, 569 
Beutler F., et al., 2014, Mon.Not.Roy.Astron.Soc., 443, 1065 
Carron J., 2013, Astron. Astrophys., 551, A88 
Carron J., Wolk M., Szapudi I., 2015, Mon.Not.Roy.Astron.Soc., 
453, 450 

Cohn J. D., 2006, New Astron., 11, 226 

Colless M., et al., 2001, Mon.Not.Roy.Astron.Soc., 328, 1039 
Crocce M., Cabre A., Gaztanaga E., 2011, 

Mon.Not.Roy.Astron.Soc., 414, 329 
Dawson K. S., et ah, 2013, Astron. J., 145, 10 
Dawson K. S., et ah, 2015, preprint (arXiv: 1508.04473) 
Dodelson S., Schneider M. D., 2013, Phys.Rev., D88, 063537 
Ellis R., et al., 2014, Publ.Astron.Soc.Jap., 66, R1 
Feldman H. A., Kaiser N., Peacock J. A., 1994, Astrophys.J., 426, 
23 

Hamilton A., 1997, Mon.Not.Roy.Astron.Soc., 289, 295 


MNRAS 000, 1-17 (2016) 


Gaussian anisotropic clustering covariance 13 


Hamilton A. J., Rimes C. D., Scoccimarro R., 2006, 

Mon.Not.Roy.Astron.Soc., 371, 1188 
Hartlap J., Simon P., Schneider P., 2006, Astron.Astrophys. 

Hill G., Gebhardt K., Komatsu E., Drory N., MacQueen P., et ah, 
2008, ASP Conf.Ser., 399, 115 
Kaiser N., 1987, Mon.Not.Roy.Astron.Soc., 227, 1 
Kalus B., Percival W., Samushia L., 2015, preprint 

(arXiv:1504.03979) 

Kaufman G., 1967, Some Bayesian Moment Formulae, Report 
No. 6710, Gentre for Operations Research and Econometrics. 
Gatholic University of Louvain, Heverlee, Belgium 
Kazin E. A., Sanchez A. G., Blanton M. R., 2012, 

Mon.Not.Roy.Astron.Soc., 419, 3223 
Labatie A., Starck J. L., Lachieze-Rey M., 2012, Astrophys. J., 
760, 97 

Laureijs R., et ah, 2011, preprint (arXiv: 1110.3193) 

Levi M., et ah, 2013, preprint (arXiv: 1308.0847) 

Lewis A., Challinor A., Lasenby A., 2000, Astrophys. J., 538, 473 
Li Y., Hu W., Takada M., 2014, Phys.Rev., D89, 083519 
Manera M., et al., 2012, Mon. Not. Roy. Astron. Soc., 428, 1036 
O’Connell R., Eisenstein D., Vargas M., Ho S., Padmanabhan N., 
2015, preprint (arXiv: 1510.01740) 

Padmanabhan N., White M. J., 2008, Phys.Rev., D77, 123540 
Paz D. J., Sanchez A. G., 2015, preprint (arXiv: 1508.03162) 
Peacock J., Smith R., 2000, Mon.Not.Roy.Astron.Soc., 318, 1144 
Pearson D. W., Samushia L., 2015, preprint (arXiv: 1509.00064) 
Percival W. J., Ross A. J., Sanchez A. G., Samushia L., Burden 
A., et al., 2014, Mon.Not.Roy.Astron.Soc., 439, 2531 
Pope A. C., Szapudi L, 2008, Mon. Not. Roy. Astron. Soc., 389, 
766 

Reid B. A., Seo H.-J., Leauthaud A., Tinker J. L., White M., 
2014, Mon. Not. Roy. Astron. Soc., 444, 476 
Ross A. J., Percival W. J., Carnero A., Zhao G.-b., Manera M., 
et al., 2013, Mon.Not.Roy.Astron.Soc., 428, 1116 
Rotenberg M., Bivins R., Metropolis N., Wooten J., 1959, The 
3-j and 6-j Symbols. Cambridge, Mass. 

Samushia L., Reid B. A., White M., Percival W. J., Cuesta A. J., 
et al., 2014, Mon. Not. Roy. Astron. Soc., 439, 3504 
Sanchez A. G., Baugh C. M., Angulo R., 2008, 

Mon.Not.Roy.Astron.Soc., 390, 1470 
Sanchez A. G., Kazin E. A., Beutler F., Chuang C.-H., Cuesta 
A. J., et al., 2013, Mon.Not.Roy.Astron.Soc., 433, 1202 
Sanchez A. G., et al., 2014, Mon. Not. Roy. Astron. Soc., 440, 
2692 

Schneider M. D., Cole S., Frenk C. S., Szapudi I., 2011, Astro- 
phys.J., 737, 11 

Scoccimarro R., Zaldarriaga M., Hui L., 1999, Astrophys.J., 527, 
1 

Scoccimarro R., Sheth R. K., Hui L., Jain B., 2001, Astrophys.J., 
546, 20 

Smith R. E., Scoccimarro R., Sheth R. K., 2008, Phys. Rev., D77, 
043525 

Springel V., 2005, Mon.Not.Roy.Astron.Soc., 364, 1105 
Springel V., White S. D., Tormen G., Kauffmann G., 2001, 
Mon.Not.Roy.Astron.Soc., 328, 726 
Sun L., Wang Q., Zhan H., 2013, Astrophys.!., 777, 75 
Takada M., Hu W., 2013, Phys.Rev., D87, 123504 
Taylor A., Joachimi B., 2014, Mon. Not. Roy. Astron. Soc., 442, 
2728 

Taylor A., Joachimi B., Kitching T., 2013, 

Mon.Not.Roy.Astron.Soc., 432, 1928 
Tinker J. L., Kravtsov A. V., Klypin A., Abazajian K., Warren 
M. S., et al., 2008, Astrophys.!., 688, 709 
White M., Reid B., Chuang C.-H., Tinker J. L., McBride C. K., 
et al., 2015, Mon.Not.Roy.Astron.Soc., 447, 234 
York D. G., et al., 2000, Astron. !., 120, 1579 
Zheng Z., Coil A. L., Zehavi I., 2007, Astrophys.!., 667, 760 


de Jong R. S., Barden S., Bellido-Tirado O., Brynnel !., Chiappini 
C., et al., 2014, Proc.SPIE Int.Soc.Opt.Eng., 9147, 91470M 
de Putter R., Wagner C., Mena O., Verde L., Percival W., 2012, 
JCAP, 1204, 019 


MNRAS 000, 1-17 (2016) 


14 J. N. Grieb et al. 




Figure Al. Left panel. Same comparison of the PS multipoles dispersion as in the lower panel of Fig. 4 but including also the predictions 
for (Tpg and crp^ for comparison. Right panel. Likewise inclusion of the predictions for and crjg to the comparison plot in the lower 
panel of Fig. 9. For brevity, only the results for the large binning scheme are shown here. 


APPENDIX A: PROOFS OF THE BINNED COVARIANCE RELATIONS 


In this section, we present the full derivation of the various formulae given in section 2.3 for the covariance of binned anisotropic 
clustering measurements. 

We assume Nh wavenumber bins centred at ki with width Afe. In the following we will average equation (14) over such 
bins, which have a volume 14^ = 47r[(fci + Ak/2)^ — {ki — AA:/2)^]/3. This is easily performed by integration in spherical 
coordinates, so that 3j:,{k — k') = Suik — k')5u{6 — d')5u{4> — (j)')/ sin0), and we will make use of the following: 
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where 14 is the volume of the galaxy sample. 


Al PS Multipoles 

Replacing P{k,fi) by its multipole expansion, P{k,fi) = Pi{k)jCe{fi), in the definition of the multipole covariance in 
equation (15), we find 
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At this point, we use the expansion of a product of two Legendre polynomials by use of Wigner 3j-symbols (Rotenberg et al. 
1959), 
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Inserting this back into the expression for o-J n gives 
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where Seo reflects that the shot noise contributes only to the monopole term. 
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First, we consider PS Legendre moments measured in the wavenumber bins defined above. Their theoretical covariance 
matrix can be obtained by the following integration over the bins: 
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(A6) 


This equation expressed in spherical coordinates as in equation (Al) reads 


Cov \Pl 




2 (2^1 + 1) (2^2 + 1) (27r) 

11,14, H 


4 pki+i^^k/l j-L 

/ k^dk 

Jki-Ak/2 J-1 


dfj, 


P{k, fi) + — 


4^1 (/r) 


where we already evaluated the Dirac delta functions. Using the definition of this yields 
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which corresponds to the result presented in equation (16). 
In particular, the monopole covariance will be given by 
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where we made use of the orthogonality of the Ci{u)- This relation could also have been derived by setting = £2 = 0 in the 
Wigner 3j-symbols in equation (16). This shows that the naive guess that monopole covariance would only be given by the 
monopole power spectrum is not correct. Instead, the monopole covariance has contributions from all higher-order multipoles. 

Note that, as implied by equation (16), the covariance of higher-order multipoles does not vanish even though the input 
power spectrum is only modelled up to £ = 4. As shown in the left panel of Fig. Al, the dispersion crp^ increases with larger 1. 


A2 PS Clustering Wedges 


For the theoretical covariance matrix of the PS clustering wedges, the integration over the wavenumber bins can be 

written as 
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by use of a similar integration as in equation (Al) restricted to the volumes of each wedge, AjiVki- The evaluation of the 
Dirac delta functions in Cov[P(fe), P(fc')] yields 
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We already presented this result in equation (17). An additional factor of 2 comes from the fact that the symmetry in /r —>■ —/t 
has been used to simplify the integration over the range into a single contiguous interval. 


A3 2PCF Multipoles 

Let us first consider the covariance of the (unbinned) 2PCF multipole moments which we will derive from the results in Fourier 
space by use of equation (5): 

Cov [ffi(s),,fr2(s')] = j d^k j dJk' Cov [Pri(fe),P(!2(fc')] e''' ® e, (A 14 ) 

where the covariance of the power spectrum multipoles is given by 

Cov[PtAk),PtAk')\ ={2nfSu{k-k')altAk)- (A15) 

The resulting expression for the covariance of 2PCF multipoles, 

Cov [ffi(s),^r 2 (s')] = J e’"'" d^k, (A 16 ) 
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Figure A2. Convergence of the dispersion of two configuration space wedges (left panel) and three wedges (right panel) with the 
multipole order included in equation (19). For two wedges, the sum can be safely truncated after £ = A without changing the results on 
a significant level. In case of three wedges, the convergence at an appropriate level is reached by including terms up to £ = 6. 


can be simplified by use of the same transformation as in equation (6), 

■£l+i2 poo 

Cov[^e^{s),^e^{s')] = ae^t^{k) ji^{ks) jt^iks ) dk. (A17) 

This equation was already given in an almost identical form by White et al. (2015), but the authors did not explicitly indicate 
the shot-noise contribution that enters the variance per power spectrum mode, o-\i^{k). 

The binned covariance of 2PCF multipole moments is then given by the volume average of the spherical Bessel functions 
over the distance bins as already presented in equation (18), 

1—^2 ^ OO 

Cov a'i^i^(k)jt^{ksi)]i^{ksj)dk, (A18) 

where we made use of the spherical Bessel functions average over the distance bin around Si with volume Vs- = 
47r - Si,min) /3 defined by 

jt(ksi) = — / ji{ks)ds. (A19) 

Vsi J si-As/2 

In contrast to the power spectrum case, the dispersion of higher-order multipole moments of the correlation function does 
not increase with 1. In fact, we show in the right panel of Fig. Al that and follow a very similar trend than the 
hexadecapole deviation. 


A4 2PCF Clustering Wedges 

We consider 2PCF clustering wedges as dehned in equation (7), and denote the measurement of in the distance bin 

around Si by 5),. Using the relation between wedges and Legendre moments given in equation (8), we can find the covariance 
of the configuration space wedges: 
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By replacing the covariance of the multipoles, we finally arrive at 

P _ I _^ j^i+^2 _ _ roo 
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(A21) 


as already presented in equation (19), where the wedge-averaged Legendre polynomials are given by equation (9). 
Evaluating the integral in that expression yields 


C. 1,11 = ^ [Ci+i{p d- Ap) - Ce-i{p + Ap) - Ce+i{p) + Ce-i{p)] ■ 


(A22) 


The sum over £i and £2 in equation (19) converges quickly. As shown in Fig. A2, the two wedge dispersion is converged 
after including only those contributions up to £i ,£2 < 4; the dispersion for three wedges needs only to also include = 6 
in the sum in order to achieve convergence at an appropriate level. 
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Figure Bl. Left panel: The values and error bars of the diagonal entries of the unbiased inverse data covariance matrix defined in 
equation (Bl), of the configuration space multipole measurements where the error haven been determined with the jackknife method. 
The dashed lines show the inverse of the theoretical covariance matrix prediction with the smoothed input power spectrum. Center 
panel: The full precision matrix of the 2PCF multipoles, Rf^^^{ki^kj) = Si) (normalized 

by the theoretical prediction), showing a complex structure and important ‘mixing terms’. Lower panel: Cut through the precision matrix 
for = 0 at Si = 112.5 Mpc. The contamination from physical effects or noise not accounted for by our modelling and is within 
the error bars which are very large due to the inversion of a noise-contaminated covariance matrix. All panels show the results for 
30 h~^ Mpc < s < 180 h~^ Mpc using As = 10 h~^ Mpc 


APPENDIX B: THE PRECISION MATRIX 


Evaluation of the Gaussian likelihood via the function given in equation (11) needs an estimate of the precision matrix. As 
pointed out in section 3.3, the covariance matrix C estimated from an ensemble of realizations is affected by noise, resulting 
in a biased estimate of the precision matrix (Hartlap et al. 2006). This bias can be removed by a rescaling of the inverse 
covariance matrix, 

where D = (Bl) 


In order not to be noise-dominated, we require D < 0.5. The only binning scheme (using three projections in jj.) that fulhls this 
requirement is the ‘large’ scheme in configuration space (As = 15 h~^ Mpc, and 12 bins per wedge/multipole measurement). 
For brevity, we only discuss the multipole results, = (1 — D) (C^)“^, where has elements Sj). Further, the 

axis-averaging of the measurements of the redshift-space two-point statistics is abandoned because equation Bl would need 
to be modified in that case. 

In analogy to the covariance notation, we write the elements of the precision matrix as , Sj). The theoretical 

prediction, Sj), is obtained from the Gaussian prediction for the multipole covariance using the smoothed input 

power spectrum (and by setting D = 0 in equation (Bl) due to the absence of sampling noise). Here we focus only on the 
non-linear prediction for illustration. The prediction from the linear input power spectrum is very similar, as already discussed 
in section 3.3. 

The measured and predicted diagonal entries of the precision matrix are shown in Fig. Bl. The error bars are determined 
with the jackknife technique (and for each jackknife estimate, the inverse covariance matrix is rescaled with a modified 
correction factor 1 — D accounting for the removed realization by Nm —^ Nm — 1). Due to the fact that the covariance matrix is 
only poorly determined, the errors on the precision matrix are very large (up to ca. 20% after the application of the correction 
factor). Within these errors, the model predictions do not show significant deviations. 

In order to show the accuracy of the modelling for the off-diagonal terms, we plot the full precision matrix normalized by 
the theoretical prediction. 


P'0 




(B2) 


in the center panel. In order to better visualize the complex structure a cut through this matrix is shown in the right panel. 
The sub-diagonal entries are negative and ‘mirror’ the diagonal entries. A similar structure is found for terms mixing the 
monopole and quadrupole as well as the quadrupole and hexadecapole entries. Further away from the block diagonals, the 
structure is noise-dominated and deviations between data and theory are largely below the noise level. 
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